Overview¶
- Freeway driving
- Introduction
- Method of mollification
- Exact solution, CLAWPACK solution and car-following model
- Night-time driving
- Introduction
- Clustering beharvior
LWR Model¶
In 1955, Lighthill, Whitham and Richards proposed the a model, now known as LWR model for traffic flow problem. This model is a macroscopic traffic flow model.
$$\rho_t+f(\rho)_x=0,$$where $$f(\rho)=\rho v(\rho).$$
Freeway Driving Model¶
From the empirical data collected from freeway, the velocity can be modeled as:
$$v(\rho)=\left\{\begin{array}{lll} 1,&\text{if }0\leq\rho<\rho_m,\\ \gamma(\frac{1}{\rho}-1),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$Hence
$$f(\rho)=\left\{\begin{array}{lll} \rho,&\text{if }0\leq\rho<\rho_m,\\ \gamma(1-\rho),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$For this problem we use $\rho_m=0.5$ and $\gamma=0.5$.
rhom = 0.5
gamma = 0.5
def flux_freeway_origin(rho):
n = rho.size
f = zeros(n)
for i in range(n):
if rho[i] >= 0. and rho[i] < rhom:
f[i] = rho[i]
elif rho[i] >= rhom and rho[i] <= 1.:
f[i] = gamma*(1.-rho[i])
else:
print "Invalid density value!"
return f
rho1 = linspace(0.,rhom-0.01,50)
rho2 = linspace(rhom,1.,51)
figure(figsize=(10,6))
plot(rho1,flux_freeway_origin(rho1),'b')
plot(rho2,flux_freeway_origin(rho2),'b')
plot([rhom,rhom],[rhom,gamma*(1.-rhom)],'--b')
ylim(0.,0.5)
xlabel('rho')
ylabel('$f$')
title('Original flux fucntion')
Method of Mollification¶
The discontinuity can result the invalidation of the technique we learned from the class, since the mehod involves the deriavtive of the flux function. One method to redeem this is to use mollifier smooth the flux function. We use the mollifier as follows:
$$\eta_\epsilon(s)=\frac{1}{\epsilon}\eta\left(\frac{s}{\epsilon}\right),$$where
$$\eta(x)=\left\{\begin{array}{lll} C\exp\left(\frac{1}{x^2-1}\right)&\text{if }-1< x < 1,\\ 0&\text{otherwise},\end{array}\right.$$Here $C$ is a constant, $C\approx2.2522836.$
The modified flux function is the convolution of the mollifier and original function:
$$f_\epsilon=\eta_\epsilon\ast f,$$so
$$f_\epsilon(\rho)=\rho+(\gamma(1-\rho)-\rho)\int_{\rho_m-\epsilon}^\rho\eta_\epsilon(s-\rho_m)ds,$$with $0<\epsilon\ll1$.
It has been proven that, if $\eta_\epsilon(x)\rightarrow\delta(x)$ as $\epsilon\rightarrow 0$, then $f_\epsilon\rightarrow f$.
In the project, we select $\epsilon=0.01$.
def eta(x,eps):
C = 2.25228362104358101049978125556
eta = where(abs(x)<eps,C*exp(eps**2/(x**2-eps**2))/eps,0.)
return eta
def flux_freeway_modified(rho,eps):
from scipy.integrate import quad
n = rho.size
f = zeros(n)
for i in range(n):
# N = 500*(int((rho[i]-rhom+eps)/eps) + 1) # mesh cell for integration
# dx = (rho[i]-rhom+eps)/N # mesh size for integration
# intl = 0.
# for j in range(N):
# intl = intl + dx*eta((j+0.5)*dx-eps,eps)
intl = quad(eta,-eps,rho[i]-rhom,args=eps)[0]
f[i] = rho[i] + (gamma*(1.-rho[i])-rho[i])*intl
return f
%%system
mkdir _plots_ini_freeway
figno = 1
fname = '_plots_ini_freeway/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
Exact Solution, Numerical Solution and Car-Following Model¶
Overview¶
Explicit Method¶
To determine the weak solution to a non-convex scalar conservation law, in this problem, we will use Olenik entropy condition:
$$ \frac{f(\rho)-f(\rho_l)}{\rho-\rho_l}\geq \frac{f(\rho_l)-f(\rho_r)}{\rho_l-\rho_r}\geq \frac{f(\rho)-f(\rho_r)}{\rho-\rho_r} $$for all $\rho$ between $\rho_l$ and $\rho_r$.
The entropy-satisfying solution to this problem can be determined form the graph of $f_\epsilon(\rho)$. This method is called convex hull.
Numerical Method¶
We use the 2nd order Lax-Wendroff method with Van Leer limiter.
Car-Following Model¶
We only focus on each car $x_n(t)$.
The density between two car is the reciprocal of the distance: $$\rho_{n+1}(t)=\frac{1}{x_{n}(t)-x_{n+1}(t)};$$
The velocity of each car only depends on the distance: $$v_{n}(t)=v(\rho_{n}(t));$$
The position is updated by: $$x_{n}(t+\Delta t)=x_{n}(t)+v_{n}(t)\Delta t.$$
Based on these assumptions, we can build up our own ODE Solver.
# ODE Solver Initialization
def initialization(N,rhol,rhor):
return linspace(0.,-(N-1.)/rhol,N)
# Main Function of ODE Solver
def ode_solver_freeway(N,tfinal,output_number,rhol,rhor):
from copy import copy
Dt = tfinal/output_number
dt = min(0.01,Dt)
x = zeros((output_number+1,N))
rho = copy(x)
x[0] = initialization(N,rhol,rhor)
rho[0][1:N+1] = -1./diff(x[0])
rho[0][0] = rhor
m = int(Dt/dt)
xs ,rhos = copy(x[0]), copy(rho[0])
for i in range(1,output_number+1):
for j in range(m):
v = flux_freeway_origin(rhos)/rhos
xs = xs + v*dt
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
v = flux_freeway_origin(rhos)/rhos
xs = xs + v*(Dt-dt*m)
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
x[i], rho[i] = copy(xs), copy(rhos)
return x, rho
N = 80
tfinal = 40.
output_number = 80
Case 1: $\rho_r<\rho_m<\rho_l$
Exact Solution¶
We construct the smallest convex hull of the set $\{(\rho,y):\rho_r<\rho<\rho_l,\ \text{and}\ y\leq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq t,\\ \rho_r \ \ \ \text{if}\ \ x > t, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$, since the shock and characteristic speeds are equal.
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.9, \rho_r=0.2.$$
PDE Solver¶
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
ODE Solver¶
figno = 1
fname = '_plots_ode_freeway1/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
Case 2: $\gamma/(1+\gamma)<\rho_l<\rho_m<\rho_r$
Exact Solution¶
We construct the smallest convex hull of the set $\{(\rho,y):\rho_l<\rho<\rho_r,\ \text{and}\ y\geq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq -\gamma t,\\ \rho_r \ \ \ \text{if}\ \ x > -\gamma t, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$.
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.4, \rho_r=0.9.$$
PDE Solver¶
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
ODE Solver¶
figno = 1
fname = '_plots_ode_freeway2/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
Case 3: $\rho_l<\rho_m<\rho_r$ and $\gamma/(1+\gamma)\geq \rho_l$
We construct the same convex hull as Case 2. There is only a shock in this solution.
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_r \ \ \ \text{if}\ \ x > st, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_r)-f_\epsilon(\rho_l)}{\rho_r-\rho_l}.$
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.2, \rho_r=0.9.$$
PDE Solver¶
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
ODE Solver¶
figno = 1
fname = '_plots_ode_freeway3/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
Clustering of Night-time Driving Model¶
Introduction¶
When driving at night on an unfamiliar mountain road, facing with the empty road will cause the speed limitation by the uncertainty of the road condition. However if there are other cars in front of the road, it will be much easier to drive faster, since their tail lights indicate how the road twists and turns.
The velocity function for night-time driving is
$$ v(\rho)=\left\{ \begin{array}{l l} v_0 \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$Hence,
$$ f(\rho)=\left\{ \begin{array}{l l} v_0\rho \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho^2 \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1\rho(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$In the car-following model, the parameters are: $$ \rho_a=0.1,\rho_b=0.3, c=10, v_0=1. $$
rhopic = linspace(0.,1.,201)
figure(figsize=(14,6))
subplot(1,2,1)
plot(rhopic,vel_night(rhopic),'b')
xlabel('rho')
ylabel('$v$')
title('Velocity')
subplot(1,2,2)
plot(rhopic,rhopic*vel_night(rhopic),'b')
xlabel('rho')
ylabel('$f$')
title('Flux')
Clustering Behavior¶
We take uniformly spaced cars as initial density on an otherwise empty road.
In LeVeque's paper:
"The lead driver sees $\rho=0$ and drives at speed $v_0$, and then the second driver initially sees $\rho=\rho_0$ and can drive faster. But as the second car accelerates, the third driver observes a drop in density and hence this driver starts to slow down. This causes the fourth driver to go faster."
It has been proven that the uniformly spaced traffic with distance $\frac{1}{\rho_0}$ with speed $v(\rho_0)$ is unstable under pertubations.
Without Perturbation¶
N = 48
tfinal = 80
output_number = 80
figure(figsize=(15,5))
rhol = 1./12
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/12$')
rhol = 1./5
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')
rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
With Perturbation¶
With random perturbations to the speeds, 2-car platoons can no longer be observed. Instead, there are random n-car platoons, which is more reasonable to model the night-time traffic.
N = 48
tfinal = 80
output_number = 80
figure(figsize=(15,5))
rhol = 1./12
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/12$')
rhol = 1./5
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')
rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
Reference¶
- LeVeque, R. J. (2001). Some traffic flow models illustrating interesting hyperbolic behavior. Minisymposium on traffic flow, Chicago.
- Wiens, J. K., Stockie, J. M., Williams, J. F. (2013). Riemann solver for a kinematic wave traffic model with discontinuous flux. Journal of Computational Physics, 242, 1-23.
- Jiang, R., Wu, Q. S. (2007). The night driving behavior in a car-following model. Physica A: Statistical Mechanics and its Applications, 375(1), 297-306.
- Dias, J. P., & Figueira, M. (2005). On the approximation of the solutions of the Riemann problem for a discontinuous conservation law. Bulletin of the Brazilian Mathematical Society, 36(1), 115-125.
- Zhang, H. M. (2003). Anisotropic property revisited––does it hold in multi-lane traffic. Transportation Research Part B: Methodological, 37(6), 561-577.